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ABSTRACT 

We present a new model of the (3 Pictoris disk-and-planet system that simulates both the plan- 
etesimal collisions and the dynamics of the resulting dust grains, allowing us to model features and 
asymmetries in both thermal and scattered light images of the disk. Our two-part model first simulates 
the collisional and dynamical evolution of the planetesimals with the Superparticle-Method Algorithm 
for Collisions in Kuiper belts (SMACK) and then simulates the dynamical evolution of the resulting 
dust grains with a standard Bulirsch-Stoer N-body integrator. Given the observed inclination and 
eccentricity of the /3 Pictoris b planet, the model neatly ties together several features of the disk: the 
central hole in the submillimeter images, the two-disk “x”-pattern seen in scattered light, and possibly 
even the clumpy gas seen by ALMA. We also hnd that most of the dust in the (3 Pictoris system is 
likely produced outside the ring at 60-100 AU. Instead of a birth ring, this disk has a “stirring ring” 
at 60-100 AU where the high-velocity collisions produced by the secular wave launched by the planet 
are concentrated. The two-disk x-pattern arises because collisions occur more frequently at the peaks 
and troughs of the secular wave. The perturbations of the disk in this region create an azimuthally 
and vertically asymmetric spatial distribution of collisions, which could yield an azimuthal clump of 
gas without invoking resonances or an additional planet. 


1. INTRODUCTION 

The (3 Pictoris system (a 21 Myr, A6 star at 19.44 pc) 
harbors a debris disk whose angular size, brightness, and 
intriguing morphology make it one of the most-obse rved 
disks in the sky. The disk was first detected by |Au- 
mann et al. ( 1984) with the Inf r ared A stronomical Satel¬ 
lite (IKAb). Smith & Terrile (1984) took the first re¬ 
solved image of the disk using a coronagraph on the du 
Pont 2.5-m telescope at the Las Campanas Observatory, 
revealing a nearly edge-on system. Eleven years after 
the first resol ved image, more deta iled images revealed 
asymmetries (Kalas & Jewitt||1995 1, including a pattern 
they referred to as a “warp" in the inner re gion of the 
disk (Burrows et al. 1995 Heap et al. 2000). In high- 


et al. (2001). Initial measurements of the planet’s orbit 
indicated th at it might be misa ligned with the second in¬ 
clined disk (Currie et al. 2011), but modeling indicated 
that, within the observational uncertainties, the observed 
planet could still be sufficientl y inclined to the ma in disk 
to be responsible for the warp (Dawson et al. 2011). More 
recent observations have been able to constrain tne orbit 
of f3 Pic b and confirm that the planet is inclined to the 
main disk ( Chauvin et al.|[2012 Nielsen et al.||2014 ). 


Submillimeter observations of the disk show a broad 


belt of planetesimals orb iting at 94 ± 8 AU (Wilner et al. 


2011 

Dent et al. 

2014 

), orbiting on the same plane as 

the extended main disi 

c of smaller dust grains observed 


in scattered light (Golimowski et al. 2006 etc.). Sub- 


resol ution scattered l i ght images of the ffisk, th e warp re- milliriieter and infrared images of the ^ Pic disk have 


pattern in the images 

7 

(Golimowski et al. 

1006 lAhmic 

et al. 112009 

Apai et al. 

11015). 

IVlouillet et al 

. (1997) and 

later Augereau et al.|(l 

' 1 11 .. u n . iD- 

(JUl 1 were able to simulate a warp 


el uding a planet with an orbit inclined to the main disk. 


Lagrange et al. (2010) later discovered a planet, j3 Pic- 
toris b, orbiting between 8 and 15 AU from the star, and 
estimated its mass a s 9 ± 3 Mjun, in agree men t with the 
model predictions of Mouillet et al. (1997) and|Augereau 
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revealed a variety of addi t ional asy mmetries in the (3 Pic- 
toris disk. |Wahhaj et al.| (|2003D and|Telesco et al.|( 2005) 
observed a bright clump of ernission in the mid-inirared 
in the SW r egion of the di sk at a projected separation 
of ~ 50 AU. Li et al. (2012) confirmed this mid-infrared 
clump at 52 AU, but noted a spatial displacement be¬ 
tween the two epochs, which they attributed to Keplerian 
motion. The nearly edge-on viewing geometry of the (3 
Pictoris debris disk complicates attempts to unravel the 
3D structure of the disk d ue to degeneracies along the 
line-of-sight. Using ALMA, Dent et al. (2014) observed 
a clump of short-lived CO gas, also in the b W region, 
but orbiting with a true orbital distance of ~ 85 AU 
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from the star, indicating an azimuthally-asym metric re¬ 
gion of enhanced collisions. Apai et al. (2015) compared 
HST/STIS observations of the disk from epochs 15 years 
apart and were unable to detect Keplerian motion from 
any point source contributing > 3% of the disk surface 
brightness at projected separations between 3'.'0 and 5"0 
(58 and 97 AU). However, this region of the disk does 
not directly probe the CO or mid-infrared clumps. 

Interpreting these patterns in the Pic disk and in 
other debris disks requires modeling collisions between 
planetesimals in addition to the dynamics of the system 
(Stark &: KuchneT] 2009 Thebault 2012 Thebault et al. 


:^01:l||Charnoz fc Tailliletll^Ol:^! INeWd et al.|!^01d||Kral 

et al. 2013| Nesvold & Kuchner 2015 Krai et al.|2015|). In 

this paper we use the term “planetesimal” to refer to the 

parent bodies in the disk larger than 1 mm in size.The 
collisional lifetime of a ~ 10 cm parent body orbiting 
in the /3 Pic disk at 10 AU is ^ 2 x 10® yr (assuming a 
low optical depth of 10“"^), while at 100 AU, the collision 
time is ~ 6 X 10^ yr. This time scale is only three time s 
the age of the system, 21 Myr (Binks & Jeffries 2014), 
and smaller bodies will have even shorter timescales tor 
fragmenting collisions, indicating that collisions are oc¬ 
curring frequently enough to influence the evolution of 
the disk. Collisions l ikely play some ro le in the observed 
ALMA asymmetry (Dent et al. 2014), highlighting the 
need for a model incorporating planetesimal collisions. 
New models of the disk also need to show how the recent 


best-fit planet orbits from Nielsen et al. (2014) influence 
the disk morphology. 

Interpreting the scattered light images calls for new 
dynamical models of the dust in the 0 Pic disk dMouille t 
et al. 1997 Augereau et al. 2001 Dawson et ^ 2011). 
W hile e xisting models can repl icate the “butterhy asym¬ 
metry” (Kalas & Jewitt 1995) or warp induced by the 
planet, they tail to reproduce the distinct two-disk x- 
pattern seen in high-resolution scatt ered light images 
(e.g.. Fig. 9 of|Golimowski et al.|2006). The dynamics of 
the dust differ from the dynamics of the planetesimals, as 
the dust is subje ct to radiation p ressure and Poynting- 


Robertson drag (Robertson 1937). The distribution of 
the dust also depends on the spatial distribution and 
collision rate of the planetesimals that produce the dust 
via collisions, as well as the destructive collisions between 
the dust grains. 

To answer this call, we present simulations of the /3 
Pictoris debris disk using our 3D collisional disk model 
SMACK, which traces the location of dust-producing col¬ 
lisions within a disk simultaneously with the evolving 
spatial distribution of larger planetesimals (> 1 mm). 
We combine these with a dynamical model of the dust 
to create the first physical model of both the planetesi¬ 
mal collisions and the dust dynamics in this system. In 
Section [2] we describe the SMACK simulations and the 
parameters used. In Sectionj^we measure the level of col¬ 
lisional relaxation within the disk to determine whether 
collisional damping is signihcant. In Section 13 we discuss 
two kinds of spiral structures created by the slightly- 
eccentric (3 Pic b, and its implications for observations 
of asymmetries in the disk. We also discuss the origin of 
the observed central hole in the distribution of mm-sized 


model of the dust in the /3 Pic disk, created by integrat¬ 
ing the orbits of the dust grains produced in SMACK 
collisions, and in Section [3 we present simulated images 
and a radial brightness profile of the dust in disk, created 
using our dust models. The model we describe in this pa¬ 
per does not include the fragmentation of grains smaller 
than 1 mm and is therefore an incomplete model of the 
small grain population, so in Section |3 we discuss the 
implications of these limitations. Finally, in Section [3 
we summarize our results and propose future work with 
this model. 

2. SMACK SIMULATIONS OF COLLIDING 
PLANETESIMALS 

Modeling the dynamical and collisional evolution of the 
parent bodies in a disk is essential for accurately calcu¬ 
lating the effects of collisional dam ping and collisional 
erosion on the morpholog y of a disk (Nesvold et al.|2013| 
Nesvold fc Kuchner|2015 ). (We do not include a full col¬ 
lisional model of the dust in this work. See Section |8] for 
a discussion of the implications of this omission.) 

Our debris disk simulator, the Superparticle-Method 
Algorithm for Collisions in Kuiper Belts (SMACK), sim¬ 
ulates the evolution of the dynamics and the spatially- 
dependent size distribution of a disk of planetesimals un¬ 
der the influence of one or more pla nets. SMACK use s 
the N-body integrator REBOUND (Rein & Liu 2012) 


but each body in the integrator represents a superpar¬ 
ticle, a cloud of planetesimals with the same position 
and trajectory, but a range of sizes from 1 mm-10 cm. 
Each superparticle is characterized by a size distribution. 
When an overlap is detected between two superparticles, 
SMACK adjusts the size distributions and trajectories 
of the colliding superparticles to represent the outcome 
of fragmenting planetesimal collisions. Data output by 
SMACK include the time-dependent 3D density map of 
the planetesimals and a 3D map of the dust-producing 
collisions throughout the simulation time. 

We simulated the evolution of the f3 Pictoris system 
with SMACK using 100,000 superparticles for a simula- 
tion time of 21 Myr, t he age of /? Pic as measured by 
Binks fc Jeffries] (2014|). We inse rted a planet o f mass 


9 Mjup, using the best-fit orbit of 


Nielsen et al. 


(2014) 


planetesimals. In Section [5T] we discuss the 3D distribu¬ 
tion of collisions in the disk. In Section we present a 


who used the Gemini/NICI and ivIageilan/iVIagAu in- 
struments to measure the position of (3 Pic b relative to 
its star with greate r accuracy than with previous obser¬ 
vatio ns alone (e.g., |Chauvin et al.||2012| |Bonnefoy et al. 
2013| . We distributed the initial superparticle orbits 
uniformly and linearly in semi-major axis, eccentricity, 
inclination, longitude of ascending node, argument of 
pericenter, and mean anomaly, and assigned their size 
distributions to be power laws with index —3.5, normal¬ 
ized such that the initial optical depth of the disk at 
100 AU was 10“^. Note that a linear distribution in 
semi-major axis will result in a radial surface density 
distribution of r~^. This radial distribution is typic ally 
for gas disk models (see |Raymond & Cossou||2014l for 
a discussion of surface density profiles), i'he initial pa¬ 
rameters of the superparticles and planet are listed in 
Table [l| We selected a superparticle size of = 10“^ ® 
AU. The finite superparticle size limits the size of the 
features that SMACK can resolve. At the location of the 
planet, the superparticle size we chose corresponds to 
an inclination and an eccentricity of Tsp/a « 0.18° and 





















































































(3 Pictoris 


3 


TABLE 1 

Initial conditions of the superparticles and planet for 
THE simulation DESCRIBED IN SECTION [5] 


Parameter 

Superparticle Range 

Planet Value 

Semi-Major Axis 

5-200 

9.1 

Eccentricity 

0.0-0.01 

0.08 

Inclination 

0.0-0.005 

0.0452 

Long, of Ascending Node 

0-27r 

0 

Argument of Periapsis 

0-271 

0.113 

Mean Anomaly 

0-277 

0 


Note. — All angular values are given in radians. 


Tgp/a « 0.003, respectively. Note that the superparti¬ 
cle encounter timescale does not represent the collisional 
timescale in the disk. The SMACK algorithm requires 
that the superparticle encounter timescale is shorter than 
the collisional timescale. It also requires a superparticle 
size small enough that the outcome of the model does not 


depend on the exact choice of superparticle size (Nesvold 
et al. 2013). The simulation ran for ^ 270 hr of wall clock 
time on the NASA Center for Climate Simulation’s Dis¬ 
cover cluster, using a hybrid OpenMP/MPI paralleliza¬ 
tion on 48 cores. 

Although we ran the simulation for 2 1 Myr, the most 


recen t estimate for the age of the star (Binks & Jeffries 


20141, the planet may have been formed or scattered onto 


its curre nt inclined orbit mor e recently than the s tar was 
formed (Dawson et al.|2011 Currie et al. 2013), so 21 
Myr does not necessarily represent the ag e of the planet- 
star sy stem in its current configuration. |Mouillet et al. 


(1997) derived a relationship between the radial extent 
of the warp, and the system age, tage^ 


emit efficiently at millimeter wavelengths. 

Our simulations reproduce the general morphology of 
previous models of the warp created by the inclined orbit 
of th e planet (e.g., Augereau et al. 2001 Dawson et al.j 


20111. At 10 Myr after the addition ot an inclined planet 
the warp has propagated out to ~ 100 AU. The shape 
of the disk as simulated with SMACK is similar to the 
results from collision less simulations (e.g.. Figure 1 in 


Dawson et al. 2011), except that our simulated image 
exhibit s a dehcit of the ma terial represented by orange 
dots in|Dawson et al. (2011). We discuss this inner clear¬ 
ing in more detail in Section |4.2[ 

Figure shows the same simulated image at 850 /xm, 
with a linear scaling and no vertical stretch. We con¬ 
volved this image with a Gaussian beam with a full-width 
half-max (FWHM) of 12 AU to simulate an ALMA ob¬ 
servation. Clearly, our simulations reproduce the basic 
morphology seen in the ALMA image, but there are some 
consequential differences between the model and obser¬ 
vations. Comparing Fi gure]^ to the cont inuum image of 
the disk with ALMA (Dent et al. 2014), we note that 
the brightness peaks correspond to the radial extent of 
the warp, which is ~ 95 AU at 10 Myr in our simulation 
(rath er than the ^ 86 AU predicted by Mouillet et al. 


1997). In the ALMA observation, the surface brightness 
peaks at ~ 60 AU on either side of the star, indicating 
that P Pictoris b reached its current inclined orbit < 10 
Myr ago, assuming the planet mass is > 9Mjup. While 
the simulated disk exhibits a small brightness asymmetry 
in this viewing orientation, our simulation is unable to 
reproduce the ~ 15% brightness asymmetry between the 
SW and NE halves of the disk seen in the ALMA image 
(see Section for a further discussion of the brightness 
asymmetry). 


= 6.31 


'mpi 

( -Tpl 

2 f 1 

\ ^age 

.M* 

UO AU7 

5.2 yr_ 


( 1 ) 


where nipi is the planet’s mass, M* is the star’s mass, 
and rpi is the planet’s orbital radius. For example, for 
a system with our parameters, in which we place the 
planet on its inclined orbit at time t = 0, the warp should 
extend out to ^ 86 AU after 10 Myr, while the maximum 
extent at 21 Myr would be ~ 107 AU. Observations of the 
disk in scattered light and submillimeter emission tend 
to place the extent of the warp in the planetesimals and 
dust between ~ 85 — 95 AU, so we synthesized images 
from our simulations after 10 Myr of evolution rather 
than 21 I^r. 

Figure [l] shows a simulated image of the /? Pictoris 
disk at lOMyr at a wavelength of 850 fj,m. To simulate 
the image, we assumed that each planetesimal within a 
superparticle emits thermally as a spherical blackbody, 
and iteratively calculated the temperature of each plan¬ 
etesimal as a function of its distance from the star. We 
summed the contribution from each superparticle in a 
given pixel, and averaged over 50 output timesteps, or 
5 X 10^ yr, to reduce the Poisson noise in the image. 
Note that at 50 AU, this is roughly two orders of mag¬ 
nitude less than the secular timescale of ^ 4.7 x 10® yr, 
ensuring that features due to secular perturbations will 
not be smeared. Separate models of the dust dynam¬ 
ics are not needed to compare to ALMA images because 
the basic SMACK models include all the grain sizes that 


3. COLLISIONAL RELAXATION 

SMACK models have the novel ability to numerically 
explore the proposed process of collisional relaxation in 
debris disks. This process, the gradual removal of free ec¬ 
centricity and free inclination from planetesimal orbits, 
has been invoked in influential models of debris disks 


with eccentric rings (e.g., Quillen 2007 Chia ng et al. 
2009|, but not investigated numerically until 


et al. 


Nesvold 

(|2013|), in which the authors demonstrated the de- 


struction, via collisions, of a spiral density wave created 
by an eccentric planet in a disk, leaving a narrow eccen¬ 
tric debris ring. 

The morphology of the /3 Pictoris disk, however, indi¬ 
cates that collisions may not have completely damped the 
free inclinations in the disk, but the collisional damping 
in P Pic has not been studied by previous models. We 
describe here what our simulations of P Pic disk show 
about this process. 

A planet on an orbit inclined to a disk will impose 
a forced inclination on the disk’s planetesimals. A plan- 
etesimal’s inclination, i, and longitude of ascending node, 
D, can be written together as a vector with components 


p = i sin D 
<7 = 1 cos fl. 


( 2 ) 

(3) 


The planetesimal’s inclination i will process about the 
forced inclination iforced induced by the planet, such that 

1 — 1/ree U ^forcedi (4) 
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Fig. 1. — Simulated SMACK image of the 0 Pictoris disk at 850 /.tm after 10 Myr. The location of the star is indicated with the white 
X. The verti cal axis has been str etched to emphasize the warp structure. The simulated disk after 10 Myr resembles the collisionless 
simulation of|Dawson et al.| l|2011|l, but in our disk model, collisions have destroyed most of the planetesimals that have completed a full 
secular oscillaition. 



Fig. 2.— Simulated SMACK image of the 0 Pictoris disk shown in Figure]^ with linear color scaling and no vertical stretch. This 
image has bee n convolved with a Gaussian beam with a FWHM of 12 AU. Comparing this simulated image with the ALMA continuum 
observation by |Dent et al.| l|2014||, we find that we reproduce the basic morphology of the ALMA image, but our brightness peaks are 
located farther out and our simulation does not reproduce the brightness asymmetry. 


where if r ep is the free or proper inc lination of the plan- 
etesimal (Murray & Dermott 1999). Inelastic collisions 
will tend to damp the tree inclinations of the superpar¬ 
ticles. 

Plotting the inclination vectors of the planetesimals in 
p-q space can help illustrate this relationship and process. 
In Figure we plot the p-q diagrams of the superpar¬ 
ticles at various times during the simulation. The black 
arrow represents the forced inclination vector imposed 
by the planet. At t = 0, the superparticles have incli¬ 
nations uniformly distributed in a small range of values, 
and are uniformly distributed in from 0 — 27r. They 
are represented on the p-q diagram as a small grey circle. 
As the system begins to evolve, the superparticles’ in¬ 
clination process about the planet’s forced inclination at 
different rates, spreading the superparticles into a ring in 
p-q space. The magnitude of the superparticles’ inclina¬ 


tion vectors oscillates between their initial inclinations 
near ~ 0 and twice the planet’s inclination, creating a 
two-component disk, with a “secondary” disk inclined to 
the main disk by twice the planet’s inclination. 

In a completely damped system, the superparticles will 
have zero free inclination, so by Equation ([4 1, their incli¬ 
nations will all equal the forced inclination of the planet, 
leaving a single disk in the plane of the planet’s or¬ 
bit rather than the two-component disk seen in 0 Pic 
and w ith collisionless simulations (e.g., Dawson et al. 


2011). On a p-q diagram, the points representing com- 
pleteiy damped superparticles would appear clustered at 
the planet’s inclination vector. We can see from the 
right panel in Figure that the superparticles are barely 
damped at all after 21 Myr. Some superparticles (rep¬ 
resenting 3% of the mass of the system) have moved in¬ 
wards towards the planet’s inclination vector due to col- 
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lisional damping of their free inclinations; others ( 1 % of 
the system by mass) have been collisionally scattered into 
higher inclinations. However, most (96% of the mass in 
the simulation) remains within the annulus defined by 
the range of their initial inclinations, indicating that the 
gravitational perturbations of the planet dominate the 
dynamical effects of collisions, and that the system is 
not yet collisionally damped. This lack of damping is 
consistent with recent measurements of the orbit of /3 
Pic b relative to the disk with Gemini/NICI and Mag- 
ellan/MagAO, which indicate that the planet’s position 
angle lies in between the positio n angles of the main disk 
and the inner warp on the skjy (Nielsen et al. 201^. 


Following the example of Uawson et al. (|2(J11 ), we 


plotted the inclination vs. semi-major axis of the su¬ 
perparticles at 10 Myr in Figure |5 Ag ain, the results 
are sim ilar to the collisionless resmts of IDawson et al.l 
(2011). The forced inclination from the planet, if, is 
independent of radial distance to the planet and simply 
equals the planet’s inclination, if = ipi. The superparti¬ 
cles orbiting at > 170 AU are still at low inclination, the 
superparticles within ^ 100 — 150 AU have just reached 
their maximum inclination of 2if for the first time, and 
the superparticles inside ^ 100 AU are oscillating be¬ 
tween the inclination of the main disk (around 0 ) and 
2if. However, in the SMACK models, we see the effects 
of collisions in the superparticles that have been colli¬ 
sionally scattered to both lower- and higher-inclination 
{2if) orbits. Again, there is little evidence of collisional 
damping, which would manifest as a clustering of super¬ 
particles at the planet’s inclination if, beginning with the 
superparticles farthest in, closest to the planet’s orbit. 

In a single collision between planetesimals, SMACK 
decreases of the kinetic energy of the planetesimals by a 
factor of 0.1 to re present the ene rgy used to fragment the 
bodies (following |Fuj iwara| 19821. Civen our previous es¬ 
timate that the collisional ntetirne of a 10 cm body at 100 
AU in the f3 Pictoris disk is 60 Myr, this indicates that 
random velocities in the disk will decrease by a factor of 
1/e after 220 Myr, consistent with the lack of discernible 
damping during our 21 Myr simulation. This estimation 
assumes that the damping rate is constant, which will 
not be the case as the damping will decrease the col¬ 
lision rate in the disk. Further simulations are needed 
to predict the timescale of the collisional damping in /3 
Pictoris and the future morphology of the disk. 

4. SPIRAL STRUCTURE 

4.1. Simulation Results 

Just as the planet’s inclination can create a warp in 
a debris disk, secular perturbations from an eccentric 
planet will perturb the orbits of the particles in the disk. 
This effect will create a spiral structure in an initially ax- 
isymmetric disk, as planetesimals at different distances 
from the planet precess at different rates. This spiral will 
spread radially outwards from the planet’s orbit. Inte¬ 
rior to the spiral, the orbits of the planetesimals are phase 
mixed, so the spiral structure only persis ts until the spi - 
ral reaches the outer edge of the disk (Wyatt 2005a). 
Co llisionless numer i cal simulations of the S Pictoris disk 
by Mouillet et al. (1997), Matthews et al. (2014| and 
Apai et al. 1 2015) predict that an eccentric (3 Pic b will 
induce a spiral in the disk via this mechanism. 


As with inclination, the eccentricity of a particle can 
be written as a vector using the particle’s longitude of 
pericenter, w. e = (e sin-07,6 cos tu). Laplace-Lagrange 
secular theory provides an expression for a planetesimal’s 
eccentricity vector e{t) perturbed by an eccentric planet. 
If we assume that the initial eccentricity of each planetes- 
imal is zero, e( 0 ) = 0, the eccentricity of a planetesimal 
with semi-major axis a is 


= 7UU- fepiCOSWpi[l 

\H/2M 

^3/2 ('^pO 


- COS At), 


03/2WPU . ^ 


(5) 


where Cpi is the planet’s eccentricity, Wpi is the planet’s 
longitude of pericenter, 63^2 63^2 Laplace coef¬ 

ficients, and api = a/upi for a > agni, where ani is the 
semi-major axis of the planet (Murray & Dermott||1999 
Wyatt et al. 1999). The precession rate A of the plan¬ 
etesimal IS given by 


^^ap;ap/^3/2(api). 


( 6 ) 


where n is the planetesimal’s mean motion, rupi is the 
mass of the planet, M, * 
ani = 1 fo r a > Upi 


is the mass of the star, and 


Murray & Dermott 


1999 


Wyatt 


The apsiaal precession period is then given 


by 2TrfA. To first order, the planetesimal’s inclination 
vector i(t), perturbed by an inclined planet, will precess 
with the same precession rate A, and the forced incli¬ 
nation is simply equal to the inclination of the planet, 
'i'pi ■ 

i (t) = {ipi cos Up; (1 -I- cos At), ipi sin Up; (1 -|- sin At)), 

(7) 

where VLr,i is the longitude of the ascending node of the 
planet (Murray & Dermott || 1999). 

To better understand the geometry of the perturbed 
disk, we used Equations [5][^ to analytically calculate the 
effects of the forced eccentricity and inclination from (3 
Pic b on the orbits after 10 Myr using the planet param¬ 
eters listed in Table and plotted the resulting orbits 
in Figure We considered the orbits of planetesimals 
on initially circular orbits with semi-ma jor axes rangin g 


from 11 AU to 155 AU. As described in Wyatt (|2005a|) 
the orbits precess at different rates, forming a spiral den- 
sity wave extending radially outward to ^ 100 AU. Inte¬ 
rior to ~ 59 AU, the planetesimals have completed more 
than one precession period and their orbits have become 
phase mixed, while exterior to ~ 100 AU, the planetesi¬ 
mals are still orbiting with very low eccentricities. 

As we described in Section the inclinations of the 
planetesimals precess around the forced eccentricity in 
a similar manner: the planetesimals interior to ^ 59 
AU have completed at least one precession period, while 
planetesimals exterior to ~ 100 AU still have very low in¬ 
clinations (Figure]^. Rather than a spiral density wave, 
however, the secular effects of the planet’s inclination in¬ 
duce a vertical displacement wave in the planetesimals. 
In Figure!^ we mark the ascending and descending nodes 
of the orbits with blue and red x’s, respectively. The 
nodes form a double-armed spiral, indicating that the 
vertical displacement wave varies azimuthally. We will 
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t = Oyr 



4«r- 


t = 2.1x10 yr 



Fig. 3. — The p-q diagrams of all of the superparticles at the beginning and end of the simulation described in Section The black 
arrow represents the forced inclination vector due to the planet. The dashed black lines illustrate the annulus in which tne inclination 
vectors of the superparticles would process in the absence of collisions. Ninety-six percent of the mass in the simulation is within the 
annulus, indicating that minimal collisional damping has occurred by 2.1 x 10^ yr. 



a(AU) 


Fig. 4. — Inclination vs. semi-major axis of the superparticles at 10 Myr. The red lines indicate the planet’s forced inclination, if, and 
2if. This figure strongly resembles the results of [Dawson et al.|(|2011|l for a > 59 AU, but at a < 59 AU in our simulations, collisions 
dominate and scatter planetesimals to high inclination. 


refer to the spiral density wave created by the planet’s 
eccentricity as the “eccentricity spiral” and the vertical 
displacement wave created by the planet’s inclination as 


the “inclination spiral”. 

Although collisions c an eventually destro y an eccen¬ 
tricity spiral in a disk (Nesvold et al. 20131, the eccen- 
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Fig. 5. — Face-on diagram of the analytically-derived orbits of planetesimals perturbed by a planet with the orbit and mass of /3 Pic b 
after 10 Myr. The black x marks the location of the star. Orbits at different semi-major axes precess at different rates, forming a spiral 
density wave where they approach adjacent orbits. Blue and red x’s mark the locations where the derived orbits cross the 2 = 0 plane, 
indicating the ascending node s and descending nodes of each orbit, respectively. They trace a double-armed spiral that intersects the spiral 
density wave (see Section |5.1[ |. 


tricity spiral induced in the /3 Pictoris disk survives to 
21 Myr. Figurej^shows a simulated image of the face-on 
(3 Pictoris disk at 10 Myr, the same simulation shown in 
Figure A spiral structure is evident, extending radi¬ 
ally outward to ~ 100 AU, in good agreement with the 
analytically-predicted spiral shown in Figure 
The outermost portion of the spiral in Figure cor¬ 
responds to the planetesimals that have completed half 
a precession period, and have reached their eccentricity 
maximum. Since the precession rate for inclinations and 
eccentricity are equal to first order, these planetesimals 
have also reached their maximum inclination, so the out¬ 
ermost spiral in the j3 Pic disk is roughly co-located with 
the maximum radial extent of the warp. 


Wilner et al. (20111 observed the j3 Pictoris disk with 
the Submillimeter Array and detected two peaks in 1.3 
mm emission along the disk plane, which they interpreted 
as a ring or belt of larger, dust-producing planetesimals 
at 94 ±8 AU, with a deficit of mm-sized grains interior to 
the belt. Our simulation results sugge st that if the warp 
in the (3 Pic disk extend s to ~ 85 AU (Golimowski et al. 


2006 Heap et al. 2000), a spiral structure created by [3 


Pic b’s eccentricity will also extend out to ^ 85 AU. The 
“ring” of planetesimals observed by Wilner et al. (20111 


may, in fact, be this spiral structure. 


Apai et al. (2015) proposed that the spiral structure 
could contribute a small b rightness asymmetry to the 


disk. Wilner et al. (2011) noted that the SW side of 


the disk appeared brighter in their SMA o bservations, 
but th e difference was below the noise level. IDent et al.l 
(|2014|) observed the /3 Pic disk at 870 ^m and found that 
the continuum emission from the SW side of the disk is 
15% brighter, on average, than from the NE side. 

To test whether the spiral structure could explain this 
asymmetry, we calculated the simulated emission from 
the NE and SW halves of the disk, observed edge-on, 
while rotating the disk about its axis, at 1 Myr intervals 
from t = 0 to t = 10 Myr. We found that the brightness 
ratio of the NE and SW regions of the disk varies with the 
orientation of the spiral at all timesteps. The maximum 
possible brightness excess of the SW disk versus the NE 
disk peaks at ^ 18% at 1 Myr, but drops to less than 
5% after the disk evolves past 4 Myr. By 10 Myr, the 
maximum possible brightness asymmetry was only ~ 2%. 
So we infer that the brightness asymmetry due to the 
spiral may contribute to the observed brightness excess 
in the SW disk, but probably cannot be solely responsible 
for it. Note that the spiral structure propagates outwards 





























Fig. 6. — Face-on simulated image of the /3 Pictoris disk at 850 /im after 10 Myr (seen edge-on in Figure The white X indicates the 
location of the star. The disk at < 60 AU is roughly an order of magnitude fainter than the disk exterior to the spiral structure because 
collisions have destroyed the planetesimals in this central clearing. 


with time and does not orbit the star, indicating that the 
brightness asymmetry due to the spiral will not move 
to the opposite side of the star, but will instead move 
radially outwards with time. 


4.2. 


Central Clearing 

The brightness of the disk shown in Figure drops 
off interior to the ring (at radii < 59 AU) by roughly 


exterior to the spiral. Wilner et al. (2011 

observed a 

deficit of mm-sized grams at radii 94 AU 

. Dent et al. 


belt with a central clearing. However, the mechanism 
for clearing these larger grains from the interior region 
of the disk is not immediately obvious. The planet will 
clear a gap around its orbit via resonance overlap, but 
this gap will only extend out to ^ 14.5 AU in 10 Myr 
even accounting f or the effects of collisions, which tend 
to widen the gap (Nesvold & Kuchner||2015 1. 

Our model suggests that a ditterent mechanism is pro¬ 
ducing the central clearing of planetesimals in the /3 Pic 
disk. The top panel of Figure shows the normalized 
radial surface brightness for the simulated disk shown in 
Figure compared with two new simulations, each run 
with 10,000 superparticles for 10 Myr. In one simula¬ 
tion, we set the eccentricity of the planet to zero. In the 


other, we kept the planet’s eccentricity set to e = 0.08, 
but we turned off the collision simulation and let the sys¬ 
tem evolve with only dynamical perturbations. Only the 
system with an eccentric planet in a disk experiencing 
c ollisions shows a br i ghtne ss deficit interior to ~ 59 AU. 


Mustill & Wyatt (2009) showed that a planet’s sec- 
ular perturbations will stir a disk of planetesimals and 
place them on intersecting orbits. Where planetesimal 
orbits cross, collisions are more frequent, and colliding 
bodies will be gradually eroded and removed from the 
system. They derived an expression for tcross, the time 
required for a planet’s secular perturbations to cause 
two neighboring planetesimal orbits to intersect. For the 
planet and star parameters in our SMACK simulation, 
tcross ~ 10 Myr at a semi-major axis of 60 AU in the 
disk. This is illustrated further in the bottom panel of 
Figure where we plot the minimum orbit intersection 
distance (MOID) versus semi-major axis for pairs of ad¬ 
jacent orbits from Figure calculated using the method 
described in |Wisniowski ^ Rickman (2013). The last 
zero of the MUiU appears to be at 59 AU. Exterior 
to 59 AU, the MOID increases, then asymptotes to a 
value of 1.1 AU, which was the separation of the orbits 
of the planetesimals in our analytic calculations. The or¬ 
bit crossings depicted by the green curve must cause the 
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Fig. 7. — Top panel: Normalized 850 /rm surface brightness vs. radius at 10 Myr for three simulations: the disk with an eccentric planet 
(shown in Figure]^, the same simulation with the planet’s eccentricity set to zero, and the eccentric planet simulation with the collisions 
turned off. An eccentric planet in a collisional disk (black curve) creates a central brightness deficit in mm-sized bodies. Bottom panel: 
Minimum orbit intersection distance (MOID) versus semi-major axis for pairs of adjacent orbits in Figure [m Orbits in the interior region 
of the disk (< 59 AU) intersect, increasing the collision rate and creating the brightness deficit see in the bmck curve in the top panel. 


collisions that create the central clearing depic ted by the 
black curve. The orbit-crossing timescale of Mustill ^ 


Wyatt (2009) indicates that the region of orbit crossing 
(and, therefore, the central clearing) will reach 94 AU at 
~ 88 Myr. 

The two methods of dis k clearing discussed in this 
section, resonance over lap (INesvold fc Kuchn er 2015) 
and secular excitation (Mustill &: W yatt 2009 ), invoke 
two different sets of initial conditions for the disk. The 


Nesvold & Kuchner (20151 model assumes that the plan- 
etesimals in the disk have some initial eccentricity dis¬ 
tribution, such that planetesimals orbiting just outside 
the resonance overlap region of the planet (the planet’s 
“cha otic zone”) will col li de fre quently and widen the gap. 
The Mustill & Wyatt (2009) mechanism assumes that 
the disk is initially cold, and depends on secular per¬ 
turbations from the planet to excite collisions between 
planetesimals. Our simulation results for the /3 Pictoris 
system suggest t hat the mechanism described by|MustiII| 


& Wyatt (2009) dominates in this system, widening the 
cleared iiiner region to much greater radial distances than 
the resonance clearing mechanism can reach in the age 


of the system. Future simulations should explore the dif¬ 
ferences between these “hot-start” and “cold-start” disk 
models, and investigate what observed disk clearing can 
tell us about the initial conditions of the disk. 

5. SIZE AND SPATIAL DISTRIBUTION OF COLLISIONS 
5.1. Spatial Distribution of Collisions 

SMACK allows us to model the 3D spatial distribution 
of collisions between parent bodies in the disk and to 
simulate spatial maps of dust production in the disk. For 
example, the left panel of Figure]^ shows a face-on map 
of the distribution of collisions in the /3 Pic disk between 
10 and 10.5 Myr. A spiral structure is evident in the 
collision rate, roughly correspondiM to the spiral seen in 
the simulated disk image in Figure!^ However, the spiral 
in the collision distribution in Figure exhibits several 
“breaks” in its azimuthal structure that are not seen in 
Figure 

These breaks can be understood by examining the 3D 
structure of the disk, specihcally the interaction between 
the eccentricity spiral (the spiral density wave induced by 
the planet’s eccentricity) and the inclination spiral (the 
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Fig. 8. — Left panel: Face-on map of superparticle collision rate between 10 and 10.5 Myr. The collision rate map contains a broken 
spiral structure corresponding roughly to the spiral density wave in Figure (the “eccentricity spiral”). Right panel: The same map, 
with the collisions in the planet’s orbital plane masked out. The mask creates a two-armed black spiral pattern (the “inclination spiral”) 
corresponding to the red and blue x’s in Figure]^ The collision rate drops where the eccentricity spiral meets the inclination spiral, creating 
the complex azimuthal structure. Two such breaks in the collision rate spiral are shown by the white arrows. 


vertical displacement wave induced by the planet’s incli¬ 
nation). In the edge-on image of the simulated disk in 
Figured the inclination spiral appears as a vertical oscil¬ 
lation m the planetesimals with radius from the star. The 
inclination spiral also appears in Figure where blue 
and red x’s mark the ascending and descending nodes 
for each orbit. The orbital nodes form a double-armed 
spiral, extending out to ~ 100 AU. As Figure shows, 
this inclination spiral intersects with the eccentricity spi¬ 
ral as several azimuthal locations. 

What happens where these two kinds of spiral inter¬ 
sect? In the right panel of Figure we plot the same 
face-on collision map shown in the left panel, but mask 
out the pixels within 0.2 AU of the plane of the planet’s 
orbit at radii > 75 AU. The right panel shows a two¬ 
armed spiral, like the pattern the ascending and descend¬ 
ing nodes illustrated in Figurej^ By comparison with the 
left panel of Figure!^ we can see that deficits in the col¬ 
lision rate occur where the inclination spirals intersect 
the eccentricity spiral in the plane of the planet’s orbit. 
Two such locations are indicated with right arrows in 
the right panel of Figure while at least two additional, 
unmarked intersection points are also visible further in. 

One reason the collision rate drops in the plane of the 
planet’s orbit is that the density drops in this plane - not 
the surface density, but the mass density. This density 
drop is shown in Figure where we plot a cut through 
the mass density of the simulated /? Pic disk at 10 Myr in 
the the a; = 0 plane. We also show a cut through the col¬ 
lision rate in the a; = 0 plane. The density and collision 
rate are enhanced at the vertical peaks of the inclination 
wave and the troughs in midplane of the initial planetes¬ 
imals disk, and minimized where the planetesimals cross 


the plane of the planet’s orbit. When the disk is viewed 
face-on (as in Figure [^, projection effects mask this ef¬ 
fect, and the spiral density wave appears continuous. 

ALMA observations show an azimuthally-asym metric 
distribution of short-lived CO gas at ^ 85 AU (Dent 
et al.|2014). The deprojected distribution of this gas ap¬ 


pears as either two clumps of gas orbiting on either side 
of the star, or a single clump with a long tail. Two ma¬ 
jor hypotheses have been suggested for the origin of the 
ALMA asymmetry: dust production from larger plan¬ 
etesimals trapped in a resonant orbit by a hypothetical 
second planet, producing the two-clump distribution, or 


2005 

Dent 

ed by 

Jack- 


the single-c lump distribution (Telesco et al. 
et al. 112014 ) . Thi s second scenario was mode, 
son el al. (2014), who found that the timescale ot tne 
observable signature of such a breakup could be as long 
as 1 Myr, but that the resulting clump would be station- 
a ry in the sys t em. 


Dent et al. (2014) concluded that the CO in the disk 
must be replenished with a steady-state production rate 
of ~ 1.4 X 10^® kg/yr. Assuming a CO ice release frac¬ 
tion of 0.1, they estimated that the mass of solid bodies 
experiencing collisions must be ^ 1.4 x 10^® kg/yr to 
maintain this CO production rate. The total mass of 
bodies involved in catastrophic collisions in our SMACK 
simulations at 10 Myr is only 1.2 x 10^® kg/yr. However, 
our simulations only track dust-producing parent bodies 
up to 10 cm in diameter. Extrapolating a simple power- 
law size distribution with index —3.5, we estimate that 
extending our initial size distribution to include bodies 
up to at least 158 km in diameter would produce the 
required 1.4 x 10^® kg/yr of colliding mass in the sim- 
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Fig. 9.— Mass density (upper panel) and collision rate (lower panel) of the simulated 0 Pic disk at 10 Myr, cut through the x = 0 plane. 
The density and collision rate are enhanced in the crests and troughs of the secular wave. 


ulated disk. This would provide enough colliding mass 
in the entire disk to produce the observed CO, but note 
that the CO in the /3 Pictoris disk is concentrated in one 
(or two) smaller subvolumes of the disk. Including still 
larger bodies may be required to produce the observed 
level of CO. Future simulations are need ed to test this 
hypothesis. Czechowski & Mann (2007) described the 
production of gas in the jJ Pic disk via high-velocity col¬ 
lisions, based on a simple 2D empirical description of the 
collision rate. We hope to improve upon this model in a 
future paper with SMACK, which can simulate both the 
3D distribution of collisions and the collisional velocities 
in the j3 Pic disk. 


5.2. Size Distributions vs. Disk Radius 

We also examined the size distribution of the parent 
bodies as they experienced collisions in SMACK. First, 
we calculated the average size distributions of superpar¬ 
ticles orbiting at radii of 50, 100, 150, and 200 AU be¬ 
tween 10.00 and 10.01 Myr, with radius bin widths of 10 
AU. We fit a power law to each distribution and found 
that the indices varied by less than 0.3 from the initial 
index of —2.5. This relatively consistent slope indicates 
that the planetesimal size distribution varies only slightly 
with radius, so the radial mass distribution of the simu¬ 
lated disk probably serves as a reasonable proxy for the 
radial optical depth distribution at submm wavelengths. 

However, these slight variations in the size distribution 
may provide information about the effects of the secular 


perturbations on the local collision rate. In Figure [TOl 
we plot each size distribution for comparison. We have 
divided the distributions by a power law with index -2.5 
to enhance the differences between the distributions and 
normalized such that the eac h d istribution equals 1 for 
the smallest size bin. Figure apparently shows that 
the spiral density wave induced by the secular pertur¬ 
bations from the planet shifts the size distribution away 
from the classical power law of index -2.5. At 10 Myr, the 
secular wave has not yet reached 200 AU, and the size 
distribution is still close to the initial power law. How¬ 
ever, at the leading edge of the wave at 150 AU, the size 
distribution becomes steeper. At 100 AU, near the peak 
density of the spiral density wave, the size distribution 
is even steeper. But after the wave has passed, at 50 
AU, the size distributions flattens back towards the -2.5 
power law, with slight deviations at the large and small 
ends. Future collisional models should investigate this 
intriguing phenomenon. 


6. INTEGRATING THE DUST ORBITS 

The SMACK models described above predict where 
dust grains are generated in the /3 Pictoris disk and their 
initial orbits. To understand the distribution of this dust 
and to model images of the disk in scattered light, we 
fed the output of SMACK into a second N-body inte- 


technique (e.g., 

Dermott et al.||1999 Liou & Zook||1999 

Wilner et al. 

2 

J02 IMoran et al. 2004 Stark & kuch- 
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Fig. 10. — Average size distributions of the superparticles at var¬ 
ious radii in the simulated disk at 10 Myr, divided by a power law 
with index -2.5 and normalized such that the number of particles 
in the smallest bin is 1 for each distribution. The size distribution 
becomes steeper as the secular wave from the planet approaches, 
then flattens back towards the initial power law. 


ner||2008 Debes et al.||2009 Kuchner & Stark 12010 1, we 
recorded the output or the ntegrator in a 3U histogram 


recorded the output of the ntegrator in a 3U histogram 
to simulate the density distribution of the dust cloud. 

We first generated a list of all the mass production 
events that took place in the simulation between 10.00 
and 10.01 Myr. Every superparticle collision found by 
SMACK yields two such mass production events, gen¬ 
erally producing different amounts of mass with differ¬ 
ent initial orbits. We recorded the semi-major axis, ec¬ 
centricity, inclination, longitude of ascending node, ar¬ 
gument of pericenter, mean anomaly, and total mass of 
dust produced for each event. The mass produced rep¬ 
resents the mass of the fragments between 1 fim and 1 
mm, distributed in a power law with an index of —0.93 
for incremental logarithmic bins (corresponding to a size 
distribution power law with index —2.8). Note that this 
mass only represents the mass of dust produced in the 
SMACK simulation during a 10,000 yr period, so the 
results we present only illustrate a relative spatial distri¬ 
bution rather than an absolute mass of dust. See Section 
8 for further discussion of this point. 

Figure [TT] illustrates the radial distribution of the most 
signihcant mass production events, binned according to 
mass production into bins containing collisions with mass 
production from 10“^® to > 10“^^, 10“^^ to 10“^® and 
> 10“^® solar masses. This figure reveals an interesting 
result. The biggest mass production events are mostly 
located in a ring roughly from 59 AU to 100 AU from the 
star, but other collisions are spread over a wider range 
of circumstellar distance, and even dominate the mass 
production at some radii. Our results indicate that only 
46% of the dust, by mass, is produced in the planetesimal 
belt (60-90 AU). In other words, attempting to apply 
the “ birth ring” approximatio n to describe the 13 Pictoris 
disk ( Strubbe fc Chiang|2006 ) would miss major sources 
of dust outside this ring" 

As a next step, we chose a subset of these mass produc¬ 
tion events to feed to the second N-body integrator. We 
selected all 960 mass production events creating > 10“^® 
solar masses of dust. We selected a similar number of 
mass production events in the range 10“^® to > 10"i® 


solar masses of dust produced, by choosing at random 
3.5% of all the events in this range. 

It is thought that the grains that dominate images of 
collision dominated disks like the (3 Pictoris disk are those 
near the blowout size, i.e. with /3 ^ 1/2 where j3 is 
the magnitude of the force on the grain from radiation 
pressure divided by the magnitude of the force on the 
grain from stellar gravity ( Strubbe fc Chian^|2006 ). We 
therefore sampled 6 different values of (3, logarithniically 
distributed, near (3 = 1/2. These values are listed in Ta¬ 
ble The grain radius corresponding to a given /3 value 
depends on assumptions about the grain’s shape, compo¬ 
sition, and optical parameters. To illustrate the range of 
possibilities, we list three different cases in Table the 
simple geometric case with solid spherica l grains with 0% 
porosi ty, and two scenarios modeled by Augereau et al.| 
(20011: 4% ice with 98% porosity and 10% ice with 95% 
porosity. For each case, we assumed a density of 3 g/cm® 
for the rock component. 

Table lists values for the corresponding grain radii 
under three possible assumptions about the grain prop¬ 
erties, as well as the fraction of the total dust surface area 
represented by each radius bin. In total, we integrated 
the orbits of 11,520 grains, assuming the simplified case 
of solid spherical grains. Since the dynamics of grains 
with /3 < 0.1 are similar, we weighted the [3 = 0.09 grains 
more to represent all grains up to a size of 1 mm. 

Our numerical integrat or uses a standard Bulirsch- 
Stoer algorithm (see e.g.. Press et al. 20071, modified 


to include terms for radiation pressure and Poynting- 


Robertson drag 

Moran et al.||200^ 

Burns et al. 

1979 Wilner et al. 

2002 

M). It Integra 

tes 

the equation ot motion 

for a dust grain I 

Robertson||1937 

), 


cPr 

df^ 


- 5 -r--- rr-l-rv, (8) 

rpO rpO Q 


where r and v are the circumstellar position and velocity 
vectors of the grain. To the right-hand side of this equa¬ 
tion, we added the gravitational force from the orbiting 
planet, though we s aw n o evidence of planet-dust inter¬ 
actions (see Section 17.11), probably because the planet is 
located so far interior to most of the dust. We chose a 
time step of one year and set the integrator to output 
the positions of the dust grains every ^ 211.56 yr (i.e., 
10.25 planet orbits). When the grains are created, their 
initial orbits conserve their birth velocities, resulting in 
high initial eccentricities, e « /3/(l — /3), for grains cre¬ 
ated in collisions bet ween two bodies on low eccentricity 
orbits . As expected (|Wyatt||2005bt |Strubbe & Chiang 
20061, we saw little Poyntmg-Robertson evolution of the 
dust grain orbits during the simulation of this collision- 
dominated disk, except for the largest values of (3. Since 
j3 Pictoris is an A-type star, we neglected the effects of 
stellar winds. 

We ran the integra tion for the grain-grain collision time 
(Wyatt et al. 1999), under the simplifying assumption 
that grains of this size are either vaporized during a 
collision or broken into daughter grains that contribute 
negligibly to the optical depth because they are quickly 
removed by radiation pressure. We approximated the 
local grain-grain collision time as a function of circum¬ 
stellar distance by calculating the local optical depth of 
the material tracked by the superparticles, and extrap- 
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Fig. 11.— Radial distribution of the mass of dust produced by collisions in the SMACK simulation, binned according to mass. The most 
massive dust production events are confined mostly to a ring between 59-100 AU, while other collisions are spread radially through the 
disk. 


TABLE 2 

The /3 values sampled for. our dust orbit integrations 




(a) 



(b) 



(c) 


/3 

s 

m 

/a 

s 

m 

fA 

s 

m 

fA 

0.09 

21 

14 

0.26 

260 

4400 

0.42 

110 

740 

0.38 

0.14 

13 

3.6 

0.057 

160 

1100 

0.045 

67 

190 

0.048 

0.23 

8.3 

0.90 

0.082 

100 

280 

0.065 

42 

47 

0.069 

0.36 

5.2 

0.23 

0.12 

65 

69 

0.093 

27 

12 

0.099 

0.57 

3.3 

0.057 

0.17 

41 

17 

0.13 

17 

3.0 

0.14 

0.90 

2.1 

0.014 

0.25 

26 

4.4 

0.19 

11 

0.74 

0.21 


Note. — We list the corresponding radii (s, in /rm), masses (m, 
in 10“® g), and the fraction of the total dust surface area represented 
by each bin (/a) for three differe nt cases: (ai Geometric optics, solid 
spherical grains, 0% p orosity, fbllAugereau e t al.l ||2001|l case for 4% 
ice, 98% porosity, (c) |Augereau et al!] (|2UUl|l case tor 10% ice, 95% 
porosity. Note that the p = U.U9 bin mcludes grains up to 1 mm in 
diameter, resulting in much larger values of /a. 


olating down 1 fim to approximate the contribution of 
smaller grains to the local optical depth. We collected 
the output coordinates into three histograms matching 
the bins used for the SMACK simulated images: one for 
the face-on image with 2 AU by 2 AU bins, and two for 
two orthogonal views of the disk edge-on, with 2 AU by 


To compare our simulated dust distributions with im¬ 
ages of the disk in scattered light, we synthesized images 
from our dust density histograms assuming that the dust 
is illuminated by 0.5 ^m light from /3 Pictoris, and scat¬ 
ters the li ght via a Henyey-G reenstein scattering phase 


function. Stark et al. (2014) found that the disk HD 


181327 was not weil-ht by a single Henyey-Greenstein 


system in which the dust is continuously replenished. We 

scattering phase function. Likewise, 

Ahmic et al. 

(2009) 

summed together the histograms representing dust from 

modeled images the (3 Pictoris disk taken by tl 

le Ad- 


the biggest mass production events with the histograms 
representing the dust from the mass production events 
in the range 10“^® to > 10“^® solar masses, weighting 
the latter histograms by 1.0/0.035 to compensate for our 
sparse sampling of these latter mass production events. 
We also applied the weightings /a listed in case (a) shown 
in Table 


7. SIMULATED SCATTERED LIGHT MORPHOLOGY 


vanced Ca mera for Surveys (ACS) on the Hubble Space 
Telescope (Golimowski et al. 2006) using a model con¬ 
sisting of a pair of thin intersecting disks, with two dif¬ 
ferent values of the Henyey-Greenstein parameter, g. In 
our physical mo del, which is not me ant to be a true in¬ 
verse model like Ahmic et al. (2009), there is no a priori 
physical distinction between the dust in one “disk” or 
another, so we used a single value for all the dust, the 
mean of the values in Ahmic et al. (2009): g = 0.743. 
Figure [T^ shows the three simulated images, showing the 
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normalized histograms of the optical depth of the dust 
from three different viewing orientations. Each image has 
been multiplied by the circumstellar distance squared to 
highlight the faint features toward the outer edge of the 
disk by revealing the underlying particle distribution. 

The face-on distribution of the dust grains exhibits a 
spiral structure, reminiscent of the spiral structure in 
the planetesimal density distribution (Figure [6). How¬ 
ever, unlike in the planetesimal distribution, which has 
a deficit of larger grains in the inner region of the disk 
(< 59 AU), there is an enhancement of dust grains in 
the inner region. This enhancement does not represent 
material falling inwards towards the star. It arises from 
the enhanced collision rate in the ;< 59 AU region of the 
disk. The planetesimals remaining in this region collide 
violently and often, due to the phase mixing of their as¬ 
cending nodes and pericenters. While our dust model 
does not simulate the 3D collisional destruction of dust 
in the disk, our simulations are supported by a spectro¬ 
scopic detecti on of silicate grains inw ards to 6 AU in the 
j3 Pic disk by |Okamoto et al.| (|2004 1. 

By comparing the bottom panel of Figure [T^ with the 
edge-on simulated image of the planetesimals in Figure 
we note that the dust grains on inclined orbits that 
make up the warp extend to greater circumstellar dis¬ 
tances than the planetesimals as t hey are pushed out¬ 
wards by radiation pressure. Figure [T^ shows that while 
a significant amount of the dust orbiting inclined to the 
main disk is within the radial extent of the planetesi¬ 
mal warp (i.e., 95 AU), there are inclined dust grains 

orbiting out to ~ 200 AU. 

Figure [goffers a more quantitative comparison of the 
models to scattered-light observations of the /3 Pictoris 
disk. It shows an edge-on vie w of the model seen in 
the bottom panel of Figure 12 (solid black curve) com¬ 
pared to the surface brightness of the disk measured by 
the ACS with arbitrary brightness scaling. The observa¬ 
tional data (dashed curve) are the power law fits t o the 
radial surfac e brightness profile listed in Table 2 of Apai| 


et al. (2015). Models and data represent the vertically- 
averaged surface brightness, and the models depict the 
side of the disk that is along the direction of the planet’s 
apocenter. The top panel of Figure [T^ shows the slope 
of the model vs. radius fr om the star , comp ared with 

(|20I5|). Different 


the power law indices fit by Apai et al. 


observational papers cited in Apai et al.|(|2Ul5|) used dif¬ 
ferent vertical heights to derive their brightness profiles, 
but their slopes agree to a level of O.I, typically, in each 
of the four regions. For our model, we averaged over the 
vertical height of ±20 AU. 

Our simulated brightness profile appears to differ sig¬ 
nificantly from the observed scattered light profiles show 
in Figure [T^ revealing the limitations of our simple dust 
model. For example, our dust model appears to be miss¬ 
ing a significant dust component. Figure shows a 
deficit in surface brightness in our model past ^ 20 AU, 
and a glance at the ACS images of the /3 Pictoris disk 
shows that our models contain too little dust in the disk 
midplane. One possible reason for these mismatches is 
that our SMACK models used the wrong initial radial 
distribution of planetesimals. For example, to conserve 
computing time, we only included planetesimals out to 
a circumstellar distance of 200 AU. See Section [H for a 
more detailed discussion of the limitations of our dust 


model. 

Our simulated surface brightness profile in Figure 13 
exhibits breaks in the the radial surface brightness power 
law, similar to but located radially inwards from the ob¬ 
served breaks. The simulated brightness profile would 
need to be scaled radially by a factor of at least 1.35 to 
fit the gross morphology of the observed profile. Given 
that the radial structure in our simulated disk propagates 
outwards in time. Equationsuggests three possible ex¬ 
planations for this discrepancy: either (a) the age of the 
planet’s orbit may be longer than 10 Myr, (b) the mass 
of /? Pic b may be larger than our simulation value of 
9 Mjup, and/or (c) the semimajor axis of /3 Pic b’s or¬ 
bit may be greater than our simulation value of 9.1 AU. 
Considering each of these scenarios individually. Equa¬ 
tion [l] requires an age of 28 Myr, a planet mass of 25 
Mjup, or a planet semimajor axis of 15 AU to scale the 
radial structure of the simulated disk by a factor of 1.35. 
These values are clearly ruled out by observations of the 
planet’s mass and orbit and the age of the system. 

A roadblock to interpretations of the radial structure 
of the P Pic disk is the contradiction between the ob¬ 
served radial extents of the vertical waves in the dust 
and pl anetesimal populations. For example, [Wilner et al.| 
(20021 measu red the radius of th e planetesimal belt to be 
~ 94 AU but Dent et al. (2014) measured the radius of 
the planetesi mal belt as ^ 50 AU, while scattered light 
observations (Heap et al. 2000 Golimowski et al. 2006) 
measured the radial extent of the warp in the dust to 
be ~ 85 AU. Future models should attempt to reconcile 
these differences to place constraints on the age of the 
planet’s current orbit. 

7.1. No Structures Orbiting With The Planet’s Mean 

Motion 

By choosing an output time step for the integrator of 
precisely 10.25 planet orbits, we can easily use our mod¬ 
els to search for disk structures that evolve on the time 
scale of t he planet’s orbit, or at harmonics of th e planet’s 
orbit (see Wilner et al. 2002 Moran et al. 2004). We con¬ 
structed four new “stroboscopic” dust density histograms 
for this purpose: one from coordinates output when the 
planet’s mean anomaly was 0, one from coordinates out¬ 
put when the planet’s mean anomaly was 0.25, and so 
on. These histograms, like four frames of a movie, can 
reveal blobs that are otherwise smeared out in other rep¬ 
resentations. 

When we constructed these histograms and differenced 
them to search for moving patterns, we found nothing 
but Poisson noise. The histograms are subject to Poisson 
noise from the finite number of particles per pixel; with 
11520 particles and typically 700 output steps per each 
stroboscopic histogram, the pixels have Poisson noise at 
the level of about l-cr = 7% per 2 AU by 2 AU pixel in 
the face-on view. So, given the limitations of our model, 
we predict that there is probably no structure rotating 
with planet’s mean motion that is stronger than about 
21 %, even at the planet’s semi-m ajor axis. _ 

This result stands in contrast to Apai et al. (20151, who 
reported a marginal detection of a ditterence between im¬ 
ages of P Pictoris taken with the Hubble Space Tele- 
yope /Space Teles cope Imaging Spectrograph in 1997 
(Heap et al. 2000) and images taken with the same in¬ 
strument m 2012. Apai et al. (2015) interpret this dif- 
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Fig. 12.— Simulated images of the /3 Pictoris disk in scattered light, multiplied by the distance to the star squared. The dust does not 
show the same deficit in the inner region of the disk as the mm-sized grains. The edge-on view in the bottom panel reveals the two-disk 
“x”-pattern seen in HST images. 


ference, seen at 3" — 6" (58-117 AU) as a possible 50% 
perturbation to the surface density. Future observations 
of this system at smaller inner working angles and future 
dynamical models that focus more on the inner region of 
the disk should resolve this issue. 

8. LIMITATIONS OF THE MODELS 

We have used our physical numerical models to explore 
and flesh out the prevailing picture of the /? Pictoris disk 
as a disk of planetesimals and dust sculpted by the sin¬ 
gle planet on an inclined orbit. But before we continue 
our discussion of how the models compare to the obser¬ 
vations, it is important to note the limitations of our 
physical models. 

One of the major limitations of this work is the use of 
two separate models: the SMACK simulations to model 
the dynamics and dust-producing collisions of the par¬ 
ent bodies > 1 mm in size, and the non-collisional dust 
dynamics model to trace the orbits of the dust grains 
under the influence of radiative forces. These two differ¬ 
ent models were applied sequentially, and simulated the 
evolution of the two populations, parent bodies and dust 
grains, with different approximations (no radiative forces 
for the parent bodies, and extremely limited collisional 


processing for the dust grains). This also prevented us 
from considering interactions between the populations. 

SMACK provides a 3D distribution of the colli¬ 
sion rate, revealing the complex azimuthally-asymmetric 
structures shown in Section |5.1[ However, in our 
SMACK simulation, we neglect planetesimals larger than 
10 cm as a source of mass in the collisional cascade. 
Thus, while we have chosen the initial conditions for our 
SMACK models to roughly match the current optical 
depth inferred from optical images, we may not have cor¬ 
rectly modeled how the absolute flux in the disk evolved 
prior to the present day. Our SMACK models also ne¬ 
glect cratering collisions, another significant source of 


dust (Theba ult & Augereau 2007 
Krivoq2010 ). Future simulations t 


Muller et al. 2010 
lat attempt to mode 
the absolute flux of the dust in j3 Pic should include the 
larger planetesimals as a source of mass as well as crater¬ 
ing collisions. Adding second-order effects such as more 
accurate grain structure and composition to future mod¬ 
els may also improve the accuracy of the simulated dust 
production in the disk. 

Our dust models assume that collisions between small 
dust grains (< 1 mm) are completely destructive, and 
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Fig. 13.— Radial surface brightness profile of our dust model (solid black line) compared with the power-law fits (dashed line) of|Apai| 
|et al.| (|2015||. The grey line shows the surface bri ghtness profile of o ur model. The upper panel shows the slope of the radial surtace 
brightness profile for our models compared with the Apai et al. ]( |2015t fits. 


calculates these collisions simply by removing grains af¬ 
ter their collisional lifetime has expired. For the smallest 
1 p,m) dust grains, the collisional products would be 
smaller than the blowout size and immediately removed 
from the system, but this assumption ignores the effects 
of fragmenting collisions of larger 100 /tm) grains, 
whose fragments may survive but whose collisions are 
not tracked by SMACK or our dust model. We there¬ 
fore cannot accurately model the size distribution of the 
dust grains. More advanced dust simulations such as a 


suitably modified collisional grooming algorithm (Stark 
&: Kuchner|2009|) (which does not currently include dust 


grain tragn i entati on), or the dust co l lision algorithms of 
Krai et al. (20131 or Vitense et al. (2014) (which treat 
dust grain iragirientation more thoroughly), could be 
used to capture the complex evolution of the dust popu¬ 
lation. 


9. DISCUSSION AND SUMMARY 

We have developed the first dynamical model for the /? 
Pictoris disk combining the colliding planetesimals, the 
dynamics of the resulting dust grains and the best-fit or¬ 
bital parameters for the planet /3 Pic b. Here, we discuss 
several features observed of the /3 Pic disk and consider 
the possible origins of these features by comparing ob¬ 
servations to our SMACK simulations and analysis. 

1. Using smack’s ability to model the dynamical ef¬ 
fects of collisions, we showed that the free incli¬ 
nations of the planetesimals are not damped sig¬ 
nificantly by collisions in the age of the system 
(Figures]^ and 1^. Our SMACK simulations repro- 


duced the warp seen in both submillimeter ( Dent 
et al.||20T4l) and scattered light observations (|Hur- 


et al.|20(i0} Golimowski et al.' 


rows et al. 


[1995 [Heap 

2006[ etc.) of the disk, confirming that this struc¬ 

ture can be induced by the planet’s inclination, in 
agreement w i th numerical sirnulation s by [Mouillet 


et al. 


et al. 


1997 1), Augereau et al. (2001), and Dawson 


2. We showed that the spiral density wave in the mm- 
sized planetesimals induced by the planet’s eccen¬ 
tricity extends out to rou ghl y the same radial ex¬ 
tent as the warp (Figures^ and [^, and could be 
i nterpreted as a ring or belt in edge-on observations 
(Wilner et al.|[20Tl] ] 


Dent et al. ][20l4l ). 


3. Our SMACK simulations also demonstrated that 
the deficit in mm-sized grains interior to this belt 


(Wilner et al. 2011 Dent et al. 


via th e mechanism described by 


20141) is created 

Iviustill fc Wyatt] 


(|2009|) : collisions are excited by the planet’s eccen¬ 
tricity interior to the spiral density wave, eroding 
and removing material from the inner region of the 
disk and creating an observable deficit (Figure]^. 


azimuthally-asymmetric dust and gas dWahhaj 

et al. 112003 

ITelesco et al.|2005[ |Li et al.|2012t 

Dent 

et al. 2014 

) via collisions, without invoking JV 


or massive collisions (Telesco et al. 

2005 

Dent 


asymmetric colfision rate afong the spiral density 
wave (Figure]^, created by the interaction between 
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the spiral density wave (induced by the planet’s ec¬ 
centricity) and the vertical displacement wave (in¬ 
duced by the planet’s inclination). Though the 
collisions among the bodies we modeled (< 10 cm 
in size) do not themselves produce enough CO to 
match the ALMA observations, collisions among 
slightly larger parent bodies in the same disk struc¬ 
ture probably could. 

5. We predict that there are no dust structures or¬ 
biting with the planet that are detectable above 
the Poisson noise of the simulation, supporting the 
suggestion that th e 50% perturbation measured by 


effects of gas drag alone probably have had a n egligible 
effect on t he du st in the (3 Pictoris disk. However, |Lyra fc| 


Kuchner (2013) demonstrated that instabilities involving 
gas drag, photoelectric heating and streaming effects can 
cause clumps and rings to form in debris disks like /3 
Pictoris. 

What is the role of large planetesimals in shaping the 
disk? Our simulations did not track bodies larger than 
10 cm in size, but the rare collisions between these mas¬ 
sive bodies could yield important contributions to the 


observable debris (e.g . 
Krai eHIlpTs] ) . 


2014 


Jackson et al. 2014 Stark et al. 


Apai et al. (2015) may have been produced by a 
recent massive collision, since such collisions were 
not included in our simulations. 


6 . Our SMACK simulations of the dust production in 
the disk revealed that only 46% of the dust is pro¬ 
duced in the planetesimal belt (60-90 AU), indicat¬ 
ing that the “bi rth ring” approximation (Strubbe 
& Chiang 2006[) fails to account for over 507o of 


What is the origin of the NE-SW brightness asymme¬ 
try? The difference in surface brightness between the 
NE and SW wings of the disk has been observed at vari¬ 
ous wavele ngths, including th e optical, infrared, and sub- 


the mass of dust produced via collision in the disk. 
Instead of a “birth ring” at this location, the /? 
Pictoris system appears to have a “stirring ring” 
where only the high-velocity planetesimal collisions 
are concentrated (Figure 13). 


millimeter (Apai et al. 2015). But perhaps the planet’s 
orbit may be more eccentric than current measurements 
indicate, or perhaps there are additional planets in the 
system creating this asymmetry. 

The 3D planetesimal and dust model we used for the 
P Pic disk could also be applied to other interesting de¬ 
bris disk systems. AU Microscopii, for example, is an¬ 
other edge-on disk with a planetesimal belt observed in 


the submillimeter (MacGregor et al. 2013) and small - 
scale structure in the dust disk (h'itzgerafd et af.||2007). 


7. Our simulated dust distribution (Figure [T^ repro¬ 
duced the x-shaped pattern seen in high-resolution 
scattered light images of the edge-on di s k by 


A combined planetesimal and dust model are needed to 
connect the submillimeter and scattered light observa¬ 
tions, by understanding how the dynamics of the plan¬ 
etesimals in the presence of a hypothetical planet affect 


et al. 

(2015 

. This nattern arises hecaiisR more rlust elers have also applied the “birth ring approximation to 

in the peaks and troughs of the secu- the AU Mic disk (Fitzgerald et al.|2007 

MacGregor et al. 

is produced 


lar wave (Figure El- It has not been repr oduced 
by previous dynamical models of the dust (|Mouil- 


let et al.||1997||Augereau et al.||2001||Uawson et al. 

201 ip. 


20131, but we have shown that the birth ring approxima- 


Though our model appears to have explained several 
salient features of the /3 Pictoris disk, many open ques¬ 
tions remain about the physics of this complicated plan¬ 
etary system. 

When and how was the planet scattered into its cur¬ 
rent orbit? Possible mechanisms suggested by Dawson 
et al. (2011) include scattering by a second planet. The 


tion can miss over half of the sources of dust in a disk 
with a planetesimal belt, indicating that observations of 
AU Mic should be revisited with a combined planetesi¬ 
mal and dust model to account for these sources. 

The HD 15115 debris disk, often referred to as the 


timescale for this scattering could be estimated by mod¬ 
eling the speed at which the secular perturbations of the 
planet propagate outwards through the disk. However, 
the observational evidence for the radial ext ent of the sec¬ 
ular p ertu rbations appears c ontradictory. |Wilner et al. 
(2002) and Dent et al. (2014) measured brightness peaks 
in the subniilfimeter observati ons, which th e y inte rpreted 


“blue needle” o r more recently, the “grey needle” (Rodi- 
gas et al.|2012 ), exhibits a brightness asymmetry despite 
its morphological symmetr y, and was recently fo und to 
have a large central cavity ( Mazoyer et al.|[2014 |. Colli- 
sional clearing may be responsible tor the central clear¬ 
ing, and the asymmetric collision effect could create a 
brightness asymmetric in the absence of an offset be¬ 
tween the disk and star. The edge-on HD 32297 also 
exhibits a brightne ss asymmetry and a central clearing 
(Currie et al. 2012), and should be studied further with 


as the planetesimal belt. But Wilner e t al. ( 2002 ) placed 
this belt at a radius of ~ 94 AU while Dent et al. (2014) 
measured it at ^ 60 AU. Scattered light observations of 
the disk indicate that the warp in the dust distribution 


extends ou t to ~ 85 AU (Heap et al. 2000 Golimowski 


et al.||20()6|. An inverse model that includes the various 


measurements of the warp and stirring belt could use 
the radial extent of these secular perturbations from the 
planet to predict when the planet was scattered onto its 
current orbit. 

What is the role of gas in s haping the dust in the disk? 
Thebault & Augereau (2005) argued that the long-term 


a collisional model. 

The asymmetric collision effect may also be observ¬ 
able in disks that are not viewed edge-on. For exam¬ 
ple, high-resolution near-IR imaging of the HD 141569A 
disk has revealed an inner spiral feat ure inclined to th e 
outer ring, and a cleared inner region ( Biller et al.|2015 ). 
A planet on an inclined and eccentric orbit could be 
responsible for such a disk morphology, and may also 
produce observable dust clumps via to asymmetric colli¬ 
sions. High-resolution submillimeter observations of HD 
141569A (for example, with ALMA) are needed to search 
for azimuthally-asymmetric gas distributions, and simu¬ 
lations like the SMACK and dust models described in 
this work could constrain the mass and orbit of a possi¬ 
ble exoplanet orbiting in the disk. 
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